Criticality of relaxation in dislocation systems 
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Relaxation processes of dislocation systems are studied by two-dimensional dynamical simulations. 
In order to capture generic features, three physically different scenarios were studied and power-law 
decays found for various physical quantities. Our main finding is that all these are the consequence 
of the underlying scaling property of the dislocation velocity distribution. Scaling is found to break 
down at some cut-off time increasing with system size. The absence of intrinsic relaxation time 
indicates that criticality is ubiquitous in all states studied. These features are reminiscent to glassy 
systems, and can be attributed to the inherent quenched disorder in the position of the slip planes. 
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When crystalline materials arc subjected to large 
enough stresses they undergo plastic, irreversible defor- 
mation caused by the motion of dislocations. As it is 
well-known, these linear lattice defects interact via long- 
range (1/r type) stress fields Q playing a crucial role in 
several complex phenomena related to plasticity, such as 
the formation of various dislocation patterns during de- 
formation and dislocation avalanches Other sys- 
tems with long-range interactions, like gravitating par- 
ticles Q or non-ncutral plasmas Q have been intensely 
studied and found to exhibit several unique properties, 
like power-law relaxation Q. Unlike dislocations these 
systems arc Hamiltonian, still one would expect similar 
level of complexities. 

Another characteristics of dislocation systems is the in- 
herent randomness in the positions of slip planes, wherein 
individual dislocations glide easily. On the other hand, 
systems with disorder, like structural and spin glasses, 
have been the focus of much interest. They were found 
to show peculiar dynamical properties 0, such as slow 
relaxation, attributed to a wide spectrum of decay times. 
The fact that dislocation systems contains quenched dis- 
order through the glide planes raises the analogy with 
spin glasses. Glassy dynamical behavior has indeed been 
observed in simulations Q and experimentally Q for dis- 
location systems, but the phenomenon is still lacking a 
systematic study. 

The interplay of long-range interaction and disorder 
leads to complex behavior even in two-dimensional (2D) 
dislocation systems [lfjl - [l3 |. The model is a strong sim- 
plification over reality, local processes 13| were neglected, 
only justified by the richness of phenomena reported 
about in this paper. Slow relaxation was observed in 
several instances, like under constant external stress, i.e., 
creep condition. In this case a dislocation system with 
single slip exhibits the well-known Andrade-typc creep 
law 1J], with the plastic strain rate 7 P i decreasing in 
time t as 7 P i(t) ~ i~ 2 / 3 until a cut-off time t\ [lij]. It was 
also suggested that t\ tends to infinity as some critical 
stress is approached from below, hinting to an analogy 
between the yielding transition and conventional phase 



transitions [15|, [l6(. In addition, it was also reported 
that single slip random 2D dislocation systems at zero 
external stress relax to an equilibrium state slowly, with 
a relaxation time increasing with system size 17] . Such 
slow relaxation processes of dislocated crystals have also 
been observed experimentally (l4. 18|. 

In this Letter we focus on the properties of relaxation 
to equilibrium below the yield stress. In all arrangements 
studied we observe power-law decay of various quanti- 
ties, a feature due to the underlying scaling of the dis- 
location velocity distribution. Scaling is always found to 
cut off at a characteristic time increasing with system 
size, indicating that critical behavior is not limited to a 
given threshold stress suggested by earlier investigations 
lot 15 1, rather, criticality is present in all states studied. 

The system considered is the simplest representation 
of a dislocated crystal, consisting of parallel straight edge 
dislocations with parallel slip planes. Thus the problem 
is simplified to 2D. By denoting the position of the ith 
dislocation by = (xi,yi), its Burgers vector by 6; = 
Si(b, 0) (sj = ±1 is the sign of the "charge" ), the equation 
of motion of a dislocation can be written as 1151 
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Here Ti n d(7") = cos(<p) cos(2(^)r _1 is the shear stress field 
generated by an individual dislocation, T oxt denotes the 
external shear stress, and N is the total number of the 
constituent dislocations. We note that the different phys- 
ical parameters are absorbed in the length-, time-, and 
stress-scales, as we measure them in p -0 ' 5 , (pMGb 2 ) -1 , 
and Gbp 05 units, respectively, where p is the dislocation 
density, M is the dislocation mobility, and G is an elastic 
constant So, in the rest of this Letter, only these 

dimcnsionlcss units are used. 

Let us highlight the most important physical features 
of this conceptually simple model. Firstly, the pair in- 
teraction is of long-range character, since it decays with 
1/r, moreover, it exhibits a complicated angular depen- 
dence with zero average. Secondly, the equation of mo- 



2 



tion ([T]) is first order, which is due to the commonly as- 
sumed overdamped nature of dislocation motion, thus 
the system is strongly dissipativc, i.e., non-Hamiltonian. 
Thirdly, the motion is constrained, since dislocations can 
only move parallel to the x axis. This has the important 
consequence that if the initial y coordinates are chosen 
randomly then this will represent quenched disorder. As 
a result, the system never completely forgets its initial 
state, that is, it does not collapse into a global energy 
minimum, rather gets trapped into a local minimum. 
Thus, the ground state of the system is frustrated leading 
to a glassy- like dynamics [|| . 

In the first part of this Letter the relaxation of random 
dislocation systems is studied. In this scenario an equal 
number of positive (sj = 1) and negative (s* = —1) sign 
dislocations are placed randomly in a square-shaped area. 
Then, with periodic boundary conditions, the equations 
of motion (Q]) are solved numerically [ItJ until no consid- 
erable dislocation movement is observed (for an example 



simulation movie see |23|). It is noted that due to the 
periodic boundary conditions, the stress field r in d is also 
modified, for details see [lij . The simulations were re- 
peated with different initial configurations 13 000, 300, 
and 118 times for the system sizes of N = 128, 512, 
and 2048, respectively During the relaxation of the sys- 
tem, the time dependent probability distribution of the 
dislocation velocities P(v, t) was determined numerically. 
(Since there is no external stress, for symmetry reasons, 
the velocity distribution of positive and negative disloca- 
tions are equal and symmetric.) According to Fig. [TJa) 
a remarkable feature of the velocity distribution is that 
it decays to zero as v~ 3 (for a theoretical explanation see 
[20|). In addition, P(v,t) tends to a Dirac-delta function 
(corresponding to the equilibrium state) and between an 
initial to and a cut-off time t± this is described with the 
scaling law 



P(v,t)=t a f(t a v). 



(2) 



The exponent a was found to be a = 0.85(2) and the 
scaling function / can be well approximated by the form 
f(x) « A/(Bx 3 + 1). Note, that the curves plotted are 
always results of averaging over the statistical ensemble. 

In order to investigate this scaling behavior in more 
detail, the mth moment (\v(t)\ m ) of the absolute velocity 
was also determined for different m values. According to 
Eq. © 



(K*)| r 



\v\ m P(v,t)dv = C m r 



(3) 



where C m is a constant. (Because of the asymptotic 
properties of the scaling function /, the integral is finite 
only for —1 < m < 2.) Figure []Jb) shows the mea- 
sured (|i>(i)| m ) curves for different m values and system 
sizes. The fitted exponents are in agreement with Eq. ([3]). 
(With the same simulation setup the scaling regime and 
the value of a has been already reported in fl7j|.) 
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FIG. 1: (color online) Dynamics during relaxation from a 
random configuration, (a) The scaled velocity distributions 
of dislocations P v given by Eq. ^ with a = 0.85 at system 
size N = 2048, and the approximated fitted scaling function 
(see text), (b) The moments (|f(i)| m ) for different exponents 
m (identified by colors) and dislocation numbers N. The 
solid, dashed, and dotted lines correspond to N = 128, 512, 
and 2048, respectively. 



As seen in Fig.fljb) the scaling region is bounded from 
both sides. It starts at a fixed to w 0.4 and lasts till a 
cut-off time t\ increasing with increasing linear system 
size L = VN. This is in perfect accordance to what 
is usually found in systems with long-range interactions. 
The actual t\ versus L relation is analyzed below. 

In the second simulation scenario, an extra disloca- 
tion with fixed position is introduced to an already re- 
laxed configuration [23]. This can be considered as the 
prototype of an external perturbation. (The equilibrium 
dislocation configuration generated by the extra disloca- 
tion, the analogue of Debye-screening, was analyzed in 
pi . 22 1.) Due to the extra force field, the system evolves 
to a new equilibrium state. As Fig.[2]shows P(v, t) obeys 
again the same scaling law given by Eq. ([2]), now with 
a = 0.34(4). Consistently, the evolution of different ve- 
locity moments obey Eq. ([3]), and thus confirms scaling. 

In the third scenario an external constant shear stress 
r cxt is turned on to the relaxed system [231 (creep ex- 
periment first studied by Miguel et al. ll|, [l5j ) . If r ex t 
is smaller than a certain yield stress, the induced mean 
plastic strain rate 7 P i(t) decreases to zero. Figure [3Ja) 
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FIG. 2: (color online) Scaling of the velocity distribution after 
adding an extra fixed dislocation to a relaxed system with 
a — 0.34 at N = 2048. The corresponding scaling function 
(solid line) is f(x) « A/(Bx :i + Cx 2 + 1). 



shows the evolution of j p \ at different external stresses 
r oxt and system sizes N. As seen 7 p ] also has a power-law 
regime with exponent slightly decreasing with increasing 
r ex t. Like for the relaxation without external stress, the 
power decay is cut off at some time t\, not shown in 
Fig. OH a) for clarity. Wc again find that t\ is increasing 
with system size, for more details see [23]. It is, however, 
practically independent from r GX t, a conclusion different 
from what was drawn by Miguel et al. [ll|, [HI . Namely, 
they found that at a critical stress level 7 P i(i) becomes 
pure power-law without a cut-off time, whereas below it 
a stress-dependent cut-off exists. According to our inves- 
tigations for low stresses the cut-off t\ oc VN is mainly 
due to finite size. 

We now turn to the velocity distribution of positive 
and negative dislocations, denoted by P + and P~, re- 
spectively. For symmetry reasons the distribution always 
fulfills P + (v,t) = P~{—v,t), so it is enough to investi- 
gate P + , denoted hereafter by P. Wc separate it into 
a symmetric and an antisymmetric part as P = P s + P a 
(note that in the absence of external stress P a = 0, hence 
P + = P~). Then in the definition of the velocity mo- 
ments [Eq. ((3])] P has to be replaced by P s , and 7 P i is 



7pi0) = (Jts SiVi / = J 



vP a (v,t)dv. 



(4) 



Figure Efb) and (c) show the scaling of P s and P a , re- 
spectively, for T GX t = 0.11 with 

P s (v,t)=t a f(t a v) and P a (v,t) = g(lfiv), (5) 

at a = 0.32(2) and j3 = 0.32(3), numerically indistin- 
guishable in this case. Note, that no power prefactor was 
found in the scaling form of P a above. From Eqs. (|4l5p 
7pl(*) = Ci~ 2/3 , with an appropriate C constant. This 
is in complete agreement with the time evolution of 7 p i 
plotted in Fig. 03a). For different r cxt values, the same 
scaling is found with slightly different exponents. Wc 
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FIG. 3: (color online) Relaxation due to a constant external 
shear stress r cxt below the yield stress, (a) The mean plas- 
tic strain rate 7 P i versus time t curves for different r cxt and 
system sizes JV. Curves for increasing size (red, green, and 
blue) correspond to N = 128, 512, and 2048, respectively. 
The plot centers on the power-decaying regime, not resolv- 
ing the cut-off times, (b) The symmetric part of the velocity 
distribution P s at different times scaled with Eq. |[5j) with 
a = 0.32 at TV = 2048 and r cxt = 0.11. The approximate 
scaling function (solid line) is f(x) ~ A/(Bx 3 + 1). (c) The 
antisymmetric part of the velocity distribution P a at differ- 
ent times scaled with Eq. © with p = 0.32 at N = 2048 and 
T" cx t = 0.11. The approximate scaling function (solid line) is 
g(x) w C/(Dx 3 +X- 1 ). 



mention that the plastic strain 7 P i(i) has the exponent 
1 — 2/3, this varies in the range 0.3 — 0.4, and is in ac- 
cordance with the well-known power 1/3 of the Andradc 
creep [14, 15|. 

To summarize the above findings, the scaling formulas 
Eq. © seem to be generally valid. The exponents vary, 
however, a feature we attribute to the difference in the 
initial conditions and external stress field, the only prop- 
erties distinguishing the scenarios we considered. Thus, 
the exponent is not a universal, inherent characteristics 
of the dislocation system determined only by the inter- 
actions and dimensionality. 

The last issue to be elucidated is the system size de- 
pendence of the cut-off time t\. In order to determine 
ti, t a (\v(t)\) is plotted with a semilogarithmic scale in 
Fig. 2] for the first simulation scenario. As seen, the 
curves obtained can be well fitted by straight lines, so 
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FIG. 4: (color online) Rescaled mean absolute velocity 
t a { | v (t ) | } as a function of scaled time t / yN at different system 
sizes. Straight line on the semilog plot shows the exponential 
component, defining a relaxation time ti for each N. Inset 
shows the obtained t± values and the ti ~ %/iV fit. 

the time evolution of (\v\) can be described with the form 
(\v{t)\) = Ct- a exp(-t/t 1 ). Thus the cut-off time h is 
a simple relaxation time. In addition, in Fig. [4] the time 
is rescaled with L = V^/V and with this the curves corre- 
sponding to different system sizes are parallel, indicating 
that t\ is proportional to the system size L. Similar be- 
havior is found for the different order velocity moments 
introduced above. These results indicate that the relax- 
ation of (|w(£)|) can be described with the rate equation 

-M>-(2 + I)<|. W |), ,6) 

where t\ cx L. This can be interpreted as follows: at 
large enough time (t 3> ti) during the relaxation of a 
dislocation system with finite size it always gets close to 
its energy minimum and enters into an exponential re- 
laxation regime. With increasing system size, however, 
the energy landscape becomes more and more complex 
resulting in a longer time to reach the exponential relax- 
ation regime. This is in agreement with the recent results 
of Laurson et al. , who found that the dynamics slow 
down dramatically as the yield stress is approached from 
above. We would also like to highlight that although the 
time exponent a varies with conditions of the relaxation, 
the size-dependence of the cut-off time t\ is in all scenar- 
ios studied in this Letter close to linear [23]. 

In summary, power-law relaxation of dislocation sys- 
tems was observed in different scenarios. This effect may 
be attributed to the quenched random positions of the 
slip axes and the complex nature of the interactions. Re- 
markably, the scaling of the time-dependent velocity dis- 
tribution goes with different exponents depending on the 
physical setup. Scaling is cut off due to the finite size, 
so the system does not possess any inherent time scale. 
The dislocation system is, therefore, found to behave like 
a critical one in all cases considered, strongly resembling 



glassy systems. 
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